function [fdata_x, fdata_y, amp1, tau1, amp2, tau2, error] = CCIV_TauFit(x_data, y_data, lam, pmask, params, maxiter, nexp)
% CCIV_TauFit - Exponential fit of rising phase of voltage trajectory
% Ultimately does double exponential curve fit; returns table of tau(2,:) and Iinj(:) from which the user
% can extract a value
% Fits directly to current-voltage relationship from current clamp.

switch (nexp)
case 2
   model=3;	%double exponential with offset
case 1
   model = 4;
otherwise
   return
end
end
FitData(:,1)=x_data';
FitData(:,2)=(y_data)';
alpha = params(1); % the delay to start of current step.
%parameters: imax     z       V0.5 
initpar = lam;
%           vrmp    A1     T1      A2     T2       inst Rin
lbound =  [-100,  -100,      0.25,   -100,     0.020,   -20];
ubound =  [ 0,     100,      100,     100,        5,    20];
order=length(initpar);

beta = 0; % beta is not used for this fit, but we need to define it anyway

[error,lam]=curve_fitting(FitData(:,1), FitData(:,2), 'levenberg','cubic', model, order, initpar,...
   pmask, lbound, ubound, alpha, beta, maxiter);

%now plot it
x = FitData(:,1); 
z = FitData(:,2);
% f=z-(lam(1).*(x-Vr)/(1+exp((x-lam(3)).*lam(2)*alpha))+lam(4).*(x-lam(5)));
fdata_x=min(x):0.1:max(x);

%calculate more values, so the fit data looks better
x1 = fdata_x(find(fdata_x < alpha));
z1 = lam(1)+0*x1;
x2 = fdata_x(find(fdata_x >= alpha));
switch (nexp)
case 2
   z2 = lam(1) + lam(2).*(1-exp(-(x2-alpha)/lam(3))) + lam(4).*(1-exp(-(x2-alpha)/lam(5)))+lam(6);
case 1
   z2 = lam(1) + lam(2).*(1-exp(-(x2-alpha)/lam(3))) + lam(4);
end

fdata_y = cat(1,z1',z2');  


switch(nexp)
case 2
   tau1 = lam(3); tau2 = lam(5);
   amp1 = lam(2); amp2 = lam(4);
   disp(sprintf('amp1: %8.3f  tau1: %8.3f  amp2: %8.3f  tau2: %8.3f', ...
      amp1, tau1, amp2, tau2))
case 1
   amp1 = lam(2); tau1 = lam(3);
   disp(sprintf('amp1: %8.3f  tau1: %8.3f', ...
      amp1, tau1))
end

% to check the fit...
%figure;
%plot(x_data, y_data, '-b', fdata_x, fdata_y, '-r');
return;